Chapter 8 HMSC analysis

8.1 Load data

load("data/data.Rdata")
load("hmsc/fit_model1_250_10.Rdata")

8.2 Variance partitioning

# Compute variance partitioning
varpart=computeVariancePartitioning(m)

varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
   group_by(variable) %>%
   summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
   tt()
tinytable_xs9ib100ukuw9za2fpnl
variable mean sd
Random: location 37.900015 25.317903
logseqdepth 56.110626 25.796874
sex 4.937460 5.612719
origin 1.051899 1.282563
# Basal tree
varpart_tree <- genome_tree

#Varpart table
varpart_table <- varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(genome=factor(genome, levels=rev(varpart_tree$tip.label))) %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location"))))

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
varpart_tree <- varpart_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
# Add phylum colors next to the tree tips
varpart_tree <- gheatmap(varpart_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
   scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
varpart_tree <- varpart_tree + new_scale_fill()

# Add variance stacked barplot
vertical_tree <-  varpart_tree +
       scale_fill_manual(values=c("#506a96","#cccccc","#be3e2b","#f6de6c"))+
        geom_fruit(
             data=varpart_table,
             geom=geom_bar,
             mapping = aes(x=value, y=genome, fill=variable, group=variable),
             pwidth = 2,
             offset = 0.05,
             width= 1,
             orientation="y",
             stat="identity")+
      labs(fill="Variable")

vertical_tree

8.3 Posterior estimates

# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree

# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    #select(genome,sp_vulgaris,area_semi,area_urban,sp_vulgarisxarea_semi,sp_vulgarisxarea_urban,season_spring,season_winter,sp_vulgarisxseason_spring,sp_vulgarisxseason_winter) %>%
    column_to_rownames(var="genome")

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)


colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
postestimates_tree <- postestimates_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
      scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

# Add posterior significant heatmap

postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.5, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
        scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
        labs(fill="Trend")

postestimates_tree +
        vexpand(.25, 1) # expand top 

8.3.1 Origin

8.3.1.1 Positively associated genomes with domestic cats

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Positive") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum, class, family)%>%
    tt()
tinytable_83tmaz02i4zlyz0iz2qb
genome phylum class order family species value
bin_CaboVerde.50 Actinomycetota Actinomycetia Actinomycetales Actinomycetaceae NA 0.953
bin_Aruba.25 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium pseudocatenulatum 0.995
bin_Aruba.2 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium adolescentis 0.994
bin_Aruba.6 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium longum 0.985
bin_Denmark.14 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium gallinarum 0.950
bin_Aruba.47 Actinomycetota Actinomycetia Mycobacteriales Mycobacteriaceae NA 0.964
bin_Aruba.57 Actinomycetota Actinomycetia Mycobacteriales Mycobacteriaceae Corynebacterium pyruviciproducens 0.931
bin_Aruba.51 Actinomycetota Coriobacteriia Coriobacteriales Atopobiaceae Parolsenella sp900544655 0.969
bin_Denmark.44 Actinomycetota Coriobacteriia Coriobacteriales Atopobiaceae UBA7748 sp900314535 0.951
bin_Aruba.14 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella stercoris 1.000
bin_Denmark.4 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella stercoris 1.000
bin_Aruba.21 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae NA 0.999
bin_Brazil.61 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella tanakaei 0.998
bin_Aruba.39 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella sp902362275 0.992
bin_Spain.84 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella sp900555555 0.992
bin_Brazil.151 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella intestinalis 0.991
bin_Denmark.127 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella sp900548365 0.991
bin_Denmark.17 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella bouchesdurhonensis 0.983
bin_Brazil.81 Bacillota Bacilli RFN20 CAG-826 NA 0.907
bin_Brazil.109 Bacillota_A Clostridia Oscillospirales Butyricicoccaceae AM07-15 sp003477405 0.936
bin_Aruba.28 Bacillota_A Clostridia Oscillospirales Oscillospiraceae NA 0.995
bin_Malaysia.9 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Dysosmobacter welbionis 0.986
bin_Malaysia.116 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Flavonifractor plautii 0.981
bin_Malaysia.34 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Lawsonibacter sp000177015 0.954
bin_Aruba.54 Bacillota_A Clostridia Tissierellales Peptoniphilaceae NA 0.959
bin_CaboVerde.35 Bacillota_A Clostridia Tissierellales Peptoniphilaceae NA 0.930
bin_Aruba.36 Bacillota_A Clostridia Tissierellales Peptoniphilaceae Peptoniphilus_C sp902363535 0.913
bin_Aruba.52 Bacillota_A Clostridia Tissierellales Peptoniphilaceae NA 0.906
bin_Malaysia.26 Bacillota_A Clostridia Oscillospirales Ruminococcaceae NA 0.976
bin_Aruba.29 Bacillota_A Clostridia Oscillospirales Ruminococcaceae Gemmiger variabilis_C 0.941
bin_Malaysia.103 Bacillota_C Negativicutes Acidaminococcales Acidaminococcaceae Acidaminococcus intestini 0.942
bin_Malaysia.8 Bacillota_C Negativicutes Veillonellales Dialisteraceae Dialister sp900543165 0.981
bin_Denmark.60 Bacillota_C Negativicutes Veillonellales Dialisteraceae Allisonella histaminiformans 0.938
bin_Aruba.64 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera sp000417505 0.970
bin_Brazil.58 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera elsdenii 0.969
bin_CaboVerde.38 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera sp000417505 0.961
bin_Malaysia.81 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera stantonii 0.943
bin_Malaysia.96 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Caecibacter sp003467125 0.919
bin_Brazil.163 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella lascolaii 0.992
bin_CaboVerde.18 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella copri 0.981
bin_Aruba.10 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella sp900544825 0.978
bin_Brazil.5 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella sp900540415 0.929
bin_Spain.21 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola coprophilus 0.929
bin_Malaysia.18 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotellamassilia sp000437675 0.924
bin_Malaysia.117 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotellamassilia sp900541335 0.919
bin_Malaysia.77 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola sp900542985 0.913
bin_Brazil.75 Campylobacterota Campylobacteria Campylobacterales Helicobacteraceae NA 0.964
bin_Malaysia.61 Campylobacterota Campylobacteria Campylobacterales Helicobacteraceae Helicobacter_C labetoulli 0.964
bin_Brazil.128 Campylobacterota Campylobacteria Campylobacterales Helicobacteraceae NA 0.963
bin_Brazil.70 Desulfobacterota Desulfovibrionia Desulfovibrionales Desulfovibrionaceae Mailhella sp003150275 0.944
bin_Brazil.146 Pseudomonadota Alphaproteobacteria RF32 CAG-239 CAG-495 sp001917125 0.969
bin_Brazil.9 Pseudomonadota Alphaproteobacteria RF32 CAG-239 CAG-495 sp000436375 0.919
bin_Aruba.15 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae NA 0.963
bin_Brazil.51 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae Sutterella wadsworthensis_A 0.963
bin_Spain.43 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae CAG-521 sp000437635 0.959
bin_Brazil.64 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae NA 0.948
bin_Brazil.92 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae CAG-521 sp000437635 0.948
bin_Brazil.15 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae NA 0.930
bin_CaboVerde.33 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum sp900543125 0.995
bin_Brazil.82 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum succiniciproducens 0.990
bin_CaboVerde.55 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum_A thomasii 0.939

8.3.1.2 Positively associated genomes feral cats

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Negative") %>%
    arrange(value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum,class,family)%>%
    tt()
tinytable_1eojf601n6x2v5g6xp6z
genome phylum class order family species value
bin_Denmark.27 Bacillota Bacilli Lactobacillales Enterococcaceae Enterococcus_B hirae 0.000
bin_Brazil.12 Bacillota Bacilli Lactobacillales Enterococcaceae Enterococcus faecalis 0.001
bin_Brazil.170 Bacillota Bacilli Lactobacillales Enterococcaceae NA 0.014
bin_CaboVerde.24 Bacillota Bacilli Lactobacillales Enterococcaceae Enterococcus_E cecorum 0.020
bin_CaboVerde.16 Bacillota Bacilli Lactobacillales Lactobacillaceae Limosilactobacillus reuteri 0.000
bin_Denmark.56 Bacillota Bacilli Lactobacillales Lactobacillaceae Levilactobacillus brevis 0.000
bin_Malaysia.72 Bacillota Bacilli Lactobacillales Lactobacillaceae Limosilactobacillus reuteri 0.000
bin_CaboVerde.60 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus agilis 0.001
bin_Malaysia.127 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus agilis 0.001
bin_Denmark.61 Bacillota Bacilli Lactobacillales Lactobacillaceae Latilactobacillus sakei 0.003
bin_Malaysia.35 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus animalis 0.003
bin_Malaysia.4 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus ruminis 0.035
bin_CaboVerde.14 Bacillota Bacilli Lactobacillales Lactobacillaceae Lactobacillus acidophilus 0.067
bin_Aruba.18 Bacillota Bacilli Lactobacillales Streptococcaceae Streptococcus lutetiensis 0.001
bin_Brazil.22 Bacillota Bacilli Lactobacillales Streptococcaceae Streptococcus pasteurianus 0.001
bin_Denmark.117 Bacillota Bacilli Lactobacillales Streptococcaceae Lactococcus garvieae 0.001
bin_Denmark.113 Bacillota Bacilli Lactobacillales Streptococcaceae Streptococcus alactolyticus 0.002
bin_Brazil.69 Bacillota_A Clostridia Oscillospirales Acutalibacteraceae Clostridium_A leptum 0.097
bin_Denmark.93 Bacillota_A Clostridia Lachnospirales Anaerotignaceae Anaerotignum sp001304995 0.074
bin_CaboVerde.3 Bacillota_A Clostridia Clostridiales Clostridiaceae Clostridium_P perfringens 0.000
bin_Denmark.42 Bacillota_A Clostridia Clostridiales Clostridiaceae Clostridium_P perfringens 0.000
bin_Brazil.136 Bacillota_A Clostridia Clostridiales Clostridiaceae Clostridium thermobutyricum 0.041
bin_Brazil.3 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900555395 0.026
bin_Brazil.89 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900550235 0.034
bin_Denmark.63 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Anaerostipes hadrus 0.040
bin_Denmark.19 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900556835 0.052
bin_Denmark.118 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Dorea_A longicatena 0.067
bin_Malaysia.108 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.081
bin_Brazil.113 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Anaerobutyricum sp900754855 0.095
bin_Denmark.50 Bacillota_A Clostridia Peptostreptococcales Peptostreptococcaceae Peptacetobacter sp900539645 0.053
bin_Brazil.107 Bacillota_A Clostridia Monoglobales_A UBA1381 CAG-41 sp900066215 0.016
bin_Brazil.159 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_B sp900541465 0.040
bin_Brazil.140 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_A sp900543175 0.046
bin_Malaysia.6 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_B sp900545035 0.088

8.3.1.3 Estimate - support plot

estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
    dplyr::select(genome,mean)

support <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
    dplyr::select(genome,support)

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
  filter(genome %in% estimate$genome) %>%
  arrange(match(genome, estimate$genome)) %>%
  dplyr::select(phylum, colors) %>%
  unique() %>%
  arrange(phylum) %>%
  dplyr::select(colors) %>%
  pull()

inner_join(estimate,support,by=join_by(genome==genome)) %>%
    mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
    mutate(support=ifelse(mean<0,1-support,support)) %>%
    left_join(genome_metadata, by = join_by(genome == genome)) %>%
    mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
    ggplot(aes(x=mean,y=support,color=phylum))+
      geom_point(alpha=0.7, shape=16, size=3)+
      scale_color_manual(values = phylum_colors) +
      geom_vline(xintercept = 0) +
      xlim(c(-0.4,0.4)) +
      labs(x = "Beta regression coefficient", y = "Posterior probability") +
      theme_minimal()+
       theme(legend.position = "none")

8.3.1.4 Correlations

#Compute the residual correlation matrix
OmegaCor = computeAssociations(m)

# Refernece tree (for sorting genomes)
genome_tree_subset <- genome_tree %>%
        keep.tip(., tip=m$spNames) 


#Co-occurrence matrix at the animal level
supportLevel = 0.95
toPlot = ((OmegaCor[[1]]$support>supportLevel)
          + (OmegaCor[[1]]$support<(1-supportLevel))>0)*OmegaCor[[1]]$mean

matrix <- toPlot %>% 
      as.data.frame() %>%
      rownames_to_column(var="genome1") %>%
      pivot_longer(!genome1, names_to = "genome2", values_to = "cor") %>%
      mutate(genome1= factor(genome1, levels=genome_tree_subset$tip.label)) %>%
      mutate(genome2= factor(genome2, levels=genome_tree_subset$tip.label)) %>%
      ggplot(aes(x = genome1, y = genome2, fill = cor)) +
            geom_tile() + 
            scale_fill_gradient2(low = "#be3e2b",
                       mid = "#f4f4f4",
                       high = "#b2b530")+
            theme_void()

htree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
vtree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
#create composite figure
grid.arrange(grobs = list(matrix,vtree),
             layout_matrix = rbind(c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1)))

8.3.1.5 Predict responses (origin)

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = c("domestic","feral")
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "origin", 
                      non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
            predict(m, Gradient = ., expected = TRUE) %>%
            do.call(rbind,.) %>%
            as.data.frame() %>%
            mutate(origin=rep(gradient,1000)) %>%
            pivot_longer(!origin,names_to = "genome", values_to = "value")
# weights:  9 (4 variable)
initial  value 101.072331 
final  value 91.392443 
converged

8.3.1.6 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)
8.3.1.6.1 Positive associated functions at element level
# Positively associated

unique_funct_db<- GIFT_db[c(2,4,5,6)] %>% 
  distinct(Code_element, .keep_all = TRUE)


element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  filter(positive_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()
tinytable_aez3gse9x22prge1df10
trait mean p1 p9 positive_support negative_support Domain Function Element
B0103 0.008668973 0.0005965939 0.016001173 0.915 0.085 Biosynthesis Nucleic acid biosynthesis UDP/UTP
D0504 0.004647085 0.0003771070 0.009768176 0.922 0.078 Degradation Amino acid degradation Methionine
D0906 0.003822519 0.0002160812 0.008201395 0.932 0.068 Degradation Antibiotic degradation Fosfomycin
D0205 0.012204272 0.0019631234 0.022408114 0.950 0.050 Degradation Polysaccharide degradation Pectin
D0208 0.009799367 0.0014936231 0.017659224 0.939 0.061 Degradation Polysaccharide degradation Mixed-Linkage glucans
element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()
tinytable_3oyen1nuof4fzv4ydsif
trait mean p1 p9 positive_support negative_support Domain Function Element
B0219 -0.004132651 -0.009242449 -2.704579e-04 0.044 0.956 Biosynthesis Amino acid biosynthesis GABA
B0214 -0.021408344 -0.038627047 -4.942490e-03 0.059 0.941 Biosynthesis Amino acid biosynthesis Glutamate
B0204 -0.016044178 -0.032154456 -1.338553e-03 0.082 0.918 Biosynthesis Amino acid biosynthesis Serine
B0302 -0.004993799 -0.010547888 -6.852718e-04 0.032 0.968 Biosynthesis Amino acid derivative biosynthesis Betaine
B0310 -0.012631605 -0.022578453 -3.174742e-03 0.035 0.965 Biosynthesis Amino acid derivative biosynthesis Tryptamine
B0303 -0.011515181 -0.020902183 -2.246292e-03 0.066 0.934 Biosynthesis Amino acid derivative biosynthesis Ectoine
B0309 -0.007441518 -0.015073476 -1.638992e-04 0.095 0.905 Biosynthesis Amino acid derivative biosynthesis Putrescine
B0804 -0.016343376 -0.029354205 -3.591691e-03 0.049 0.951 Biosynthesis Aromatic compound biosynthesis Dipicolinate
B0603 -0.016580393 -0.031741463 -2.636490e-03 0.060 0.940 Biosynthesis Organic anion biosynthesis Citrate
B0601 -0.009240453 -0.017544737 -1.391204e-03 0.069 0.931 Biosynthesis Organic anion biosynthesis Succinate
B0401 -0.011102735 -0.022069357 -2.102256e-04 0.094 0.906 Biosynthesis SCFA biosynthesis Acetate
B0709 -0.002073624 -0.003638992 -6.474296e-04 0.033 0.967 Biosynthesis Vitamin biosynthesis Tocopherol/tocotorienol (E)
D0517 -0.004495743 -0.007874479 -1.353705e-03 0.023 0.977 Degradation Amino acid degradation Ornithine
D0508 -0.003245058 -0.007138266 -1.816237e-04 0.079 0.921 Degradation Amino acid degradation Lysine
D0512 -0.018089246 -0.035972853 -1.186182e-03 0.086 0.914 Degradation Amino acid degradation Histidine
D0903 -0.004117914 -0.009256580 -2.622164e-04 0.044 0.956 Degradation Antibiotic degradation Cephalosporin
D0908 -0.015447072 -0.027438899 -3.764076e-03 0.056 0.944 Degradation Antibiotic degradation Macrolide
D0601 -0.009517923 -0.017285585 -2.651368e-03 0.024 0.976 Degradation Nitrogen compound degradation Nitrate
D0603 -0.001981904 -0.003844551 -3.973787e-04 0.036 0.964 Degradation Nitrogen compound degradation Urate
D0610 -0.002955038 -0.004863991 -1.021881e-03 0.039 0.961 Degradation Nitrogen compound degradation Methylamine
D0611 -0.004117914 -0.009256580 -2.622164e-04 0.044 0.956 Degradation Nitrogen compound degradation Phenylethylamine
D0606 -0.005881451 -0.011632320 -9.896364e-04 0.063 0.937 Degradation Nitrogen compound degradation Allantoin
D0612 -0.001574312 -0.002955348 -1.518534e-04 0.076 0.924 Degradation Nitrogen compound degradation Hypotaurine
D0801 -0.001305868 -0.002159595 -9.074796e-05 0.008 0.992 Degradation Xenobiotic degradation Toluene
D0802 -0.001305868 -0.002159595 -9.074796e-05 0.008 0.992 Degradation Xenobiotic degradation Xylene
D0807 -0.004098974 -0.008643808 -7.646272e-04 0.043 0.957 Degradation Xenobiotic degradation Catechol
D0817 -0.004953919 -0.010591016 -6.720724e-04 0.049 0.951 Degradation Xenobiotic degradation Trans-cinnamate
D0816 -0.005871819 -0.011500014 -4.581529e-04 0.083 0.917 Degradation Xenobiotic degradation Phenylacetate
8.3.1.6.2 Negatively associated functions at element level
element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()
tinytable_jjjw1cnk4hd23g99arcg
trait mean p1 p9 positive_support negative_support Domain Function Element
B0219 -0.004132651 -0.009242449 -2.704579e-04 0.044 0.956 Biosynthesis Amino acid biosynthesis GABA
B0214 -0.021408344 -0.038627047 -4.942490e-03 0.059 0.941 Biosynthesis Amino acid biosynthesis Glutamate
B0204 -0.016044178 -0.032154456 -1.338553e-03 0.082 0.918 Biosynthesis Amino acid biosynthesis Serine
B0302 -0.004993799 -0.010547888 -6.852718e-04 0.032 0.968 Biosynthesis Amino acid derivative biosynthesis Betaine
B0310 -0.012631605 -0.022578453 -3.174742e-03 0.035 0.965 Biosynthesis Amino acid derivative biosynthesis Tryptamine
B0303 -0.011515181 -0.020902183 -2.246292e-03 0.066 0.934 Biosynthesis Amino acid derivative biosynthesis Ectoine
B0309 -0.007441518 -0.015073476 -1.638992e-04 0.095 0.905 Biosynthesis Amino acid derivative biosynthesis Putrescine
B0804 -0.016343376 -0.029354205 -3.591691e-03 0.049 0.951 Biosynthesis Aromatic compound biosynthesis Dipicolinate
B0603 -0.016580393 -0.031741463 -2.636490e-03 0.060 0.940 Biosynthesis Organic anion biosynthesis Citrate
B0601 -0.009240453 -0.017544737 -1.391204e-03 0.069 0.931 Biosynthesis Organic anion biosynthesis Succinate
B0401 -0.011102735 -0.022069357 -2.102256e-04 0.094 0.906 Biosynthesis SCFA biosynthesis Acetate
B0709 -0.002073624 -0.003638992 -6.474296e-04 0.033 0.967 Biosynthesis Vitamin biosynthesis Tocopherol/tocotorienol (E)
D0517 -0.004495743 -0.007874479 -1.353705e-03 0.023 0.977 Degradation Amino acid degradation Ornithine
D0508 -0.003245058 -0.007138266 -1.816237e-04 0.079 0.921 Degradation Amino acid degradation Lysine
D0512 -0.018089246 -0.035972853 -1.186182e-03 0.086 0.914 Degradation Amino acid degradation Histidine
D0903 -0.004117914 -0.009256580 -2.622164e-04 0.044 0.956 Degradation Antibiotic degradation Cephalosporin
D0908 -0.015447072 -0.027438899 -3.764076e-03 0.056 0.944 Degradation Antibiotic degradation Macrolide
D0601 -0.009517923 -0.017285585 -2.651368e-03 0.024 0.976 Degradation Nitrogen compound degradation Nitrate
D0603 -0.001981904 -0.003844551 -3.973787e-04 0.036 0.964 Degradation Nitrogen compound degradation Urate
D0610 -0.002955038 -0.004863991 -1.021881e-03 0.039 0.961 Degradation Nitrogen compound degradation Methylamine
D0611 -0.004117914 -0.009256580 -2.622164e-04 0.044 0.956 Degradation Nitrogen compound degradation Phenylethylamine
D0606 -0.005881451 -0.011632320 -9.896364e-04 0.063 0.937 Degradation Nitrogen compound degradation Allantoin
D0612 -0.001574312 -0.002955348 -1.518534e-04 0.076 0.924 Degradation Nitrogen compound degradation Hypotaurine
D0801 -0.001305868 -0.002159595 -9.074796e-05 0.008 0.992 Degradation Xenobiotic degradation Toluene
D0802 -0.001305868 -0.002159595 -9.074796e-05 0.008 0.992 Degradation Xenobiotic degradation Xylene
D0807 -0.004098974 -0.008643808 -7.646272e-04 0.043 0.957 Degradation Xenobiotic degradation Catechol
D0817 -0.004953919 -0.010591016 -6.720724e-04 0.049 0.951 Degradation Xenobiotic degradation Trans-cinnamate
D0816 -0.005871819 -0.011500014 -4.581529e-04 0.083 0.917 Degradation Xenobiotic degradation Phenylacetate
positive <- element_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- element_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

table <- bind_rows(positive,negative) %>%
  left_join(unique_funct_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) 

table %>%
  mutate(Element=factor(Element,levels=c(table$Element))) %>%
  ggplot(aes(x=mean, y=fct_rev(Element), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.3.1.7 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })
#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
function_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()
tinytable_tk7jnsi0uenfjp88zy7f
trait mean p1 p9 positive_support negative_support
D02 8.252038e-03 -0.0029397646 0.0196760509 0.822 0.178
B08 7.109220e-03 -0.0040068882 0.0171378858 0.770 0.230
B01 1.216320e-03 -0.0059241188 0.0081540077 0.610 0.390
S01 1.014944e-03 -0.0109145072 0.0133945495 0.551 0.449
B09 8.119340e-07 -0.0005326007 0.0004986745 0.367 0.633
# Negatively associated
function_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()
tinytable_qa560v83hjg4lsvwdyzf
trait mean p1 p9 positive_support negative_support
D08 -1.093511e-03 -0.002214746 -2.281340e-04 0.032 0.968
B03 -1.057672e-02 -0.017631148 -3.229192e-03 0.045 0.955
D06 -3.229620e-03 -0.006962169 -5.139396e-05 0.098 0.902
B04 -7.665435e-03 -0.017979073 1.779859e-03 0.149 0.851
B06 -6.778849e-03 -0.016584138 2.083246e-03 0.171 0.829
D07 -1.166506e-02 -0.029621342 4.098672e-03 0.194 0.806
D03 -4.540739e-03 -0.012618543 2.717105e-03 0.207 0.793
D05 -1.909669e-03 -0.007386594 2.833109e-03 0.221 0.779
B02 -3.789939e-03 -0.012493886 4.029648e-03 0.247 0.753
S03 -9.214706e-03 -0.031916056 1.659255e-02 0.249 0.751
D09 -1.644702e-03 -0.007831696 4.607589e-03 0.318 0.682
S02 -4.283466e-03 -0.013942829 2.781315e-03 0.319 0.681
B07 -3.643201e-03 -0.015378159 8.373293e-03 0.322 0.678
D01 -3.039883e-04 -0.004674591 4.160412e-03 0.417 0.583
B10 -1.661264e-05 -0.000312485 2.497855e-04 0.449 0.551
positive <- function_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- function_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.02,0.02)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.3.2 Sex

8.3.2.1 Negatively associated genomes with male cats

# Compute variance partitioning
varpart=computeVariancePartitioning(m)

varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
   group_by(variable) %>%
   summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
   tt()
tinytable_vmb9rx0xx39470t262i5
variable mean sd
Random: location 37.900015 25.317903
logseqdepth 56.110626 25.796874
sex 4.937460 5.612719
origin 1.051899 1.282563
# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree
# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    column_to_rownames(var="genome")

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="sexMale", trend=="Negative") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum, class, family)%>%
    tt()
tinytable_rq04k0amne019wihr37m
genome phylum class order family species value
bin_Denmark.4 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella stercoris 0.090
bin_Aruba.14 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella stercoris 0.089
bin_Brazil.91 Bacillota Bacilli RF39 UBA660 CAG-988 sp003149915 0.082
bin_Brazil.76 Bacillota Bacilli RF39 UBA660 CAG-877 sp900554305 0.017
bin_Malaysia.17 Bacillota Bacilli RF39 UBA660 NA 0.016
bin_Brazil.45 Bacillota Bacilli RF39 UBA660 CAG-877 sp900554315 0.010
bin_Malaysia.16 Bacillota_A Clostridia Oscillospirales Acutalibacteraceae NA 0.085
bin_Malaysia.78 Bacillota_A Clostridia Oscillospirales Acutalibacteraceae Ruminococcus_E bromii_B 0.081
bin_Brazil.69 Bacillota_A Clostridia Oscillospirales Acutalibacteraceae Clostridium_A leptum 0.020
bin_Brazil.83 Bacillota_A Clostridia Lachnospirales Anaerotignaceae NA 0.069
bin_Denmark.63 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Anaerostipes hadrus 0.091
bin_Spain.11 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.084
bin_Malaysia.3 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.077
bin_Malaysia.45 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Blautia_A wexlerae 0.077
bin_Brazil.1 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.071
bin_Malaysia.7 Bacillota_A Clostridia Lachnospirales Lachnospiraceae CAG-317 sp900543415 0.070
bin_Denmark.118 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Dorea_A longicatena 0.066
bin_Brazil.3 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900555395 0.065
bin_Malaysia.52 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Fusicatenibacter saccharivorans 0.065
bin_Brazil.89 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900550235 0.062
bin_Denmark.34 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Fusicatenibacter saccharivorans 0.062
bin_Malaysia.73 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.062
bin_Brazil.166 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.057
bin_Brazil.165 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Blautia_A caecimuris 0.055
bin_Denmark.19 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900556835 0.053
bin_Spain.53 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.053
bin_Brazil.54 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Dorea_B phocaeensis 0.051
bin_Brazil.157 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.046
bin_Brazil.32 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Mediterraneibacter torques 0.042
bin_Brazil.8 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Blautia_A sp900547615 0.042
bin_Spain.6 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Schaedlerella glycyrrhizinilytica 0.039
bin_Denmark.109 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Eisenbergiella sp900550285 0.034
bin_Brazil.113 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Anaerobutyricum sp900754855 0.031
bin_Brazil.17 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.030
bin_Spain.37 Bacillota_A Clostridia Lachnospirales Lachnospiraceae UMGS1472 sp900552095 0.029
bin_Denmark.52 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Ruminococcus_B gnavus 0.023
bin_Brazil.56 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Ruminococcus_A sp000432335 0.022
bin_Spain.67 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Ruminococcus_B sp900544395 0.021
bin_Brazil.28 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.019
bin_Brazil.93 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.019
bin_Malaysia.110 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Catenibacillus sp900557175 0.018
bin_Aruba.43 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Roseburia sp900548205 0.017
bin_Malaysia.108 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.017
bin_Brazil.116 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.016
bin_Malaysia.86 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.015
bin_Brazil.99 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Lachnospira sp900552795 0.013
bin_CaboVerde.61 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Roseburia sp900548205 0.013
bin_Brazil.97 Bacillota_A Clostridia Lachnospirales Lachnospiraceae UBA9502 sp900538475 0.010
bin_Malaysia.46 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Blautia sp003287895 0.010
bin_Brazil.50 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Lachnospira sp003451515 0.009
bin_Denmark.66 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Blautia sp900120295 0.009
bin_Spain.60 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.009
bin_Malaysia.125 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Blautia sp900539145 0.008
bin_CaboVerde.34 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.005
bin_Aruba.13 Bacillota_A Clostridia Lachnospirales Lachnospiraceae CAG-81 sp000435795 0.002
bin_Brazil.63 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.002
bin_Denmark.62 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Clostridium_Q sp000435655 0.002
bin_Brazil.105 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.001
bin_Denmark.20 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Enterocloster sp001517625 0.001
bin_Malaysia.98 Bacillota_A Clostridia Oscillospirales Oscillospiraceae NA 0.066
bin_Malaysia.9 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Dysosmobacter welbionis 0.065
bin_Aruba.28 Bacillota_A Clostridia Oscillospirales Oscillospiraceae NA 0.044
bin_Brazil.177 Bacillota_A Clostridia Oscillospirales Oscillospiraceae NA 0.044
bin_Malaysia.116 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Flavonifractor plautii 0.043
bin_Malaysia.34 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Lawsonibacter sp000177015 0.040
bin_Brazil.137 Bacillota_A Clostridia Oscillospirales Oscillospiraceae NA 0.027
bin_Aruba.52 Bacillota_A Clostridia Tissierellales Peptoniphilaceae NA 0.094
bin_Aruba.31 Bacillota_A Clostridia Oscillospirales Ruminococcaceae Negativibacillus sp000435195 0.080
bin_Malaysia.102 Bacillota_A Clostridia Oscillospirales Ruminococcaceae NA 0.080
bin_Malaysia.30 Bacillota_A Clostridia Oscillospirales Ruminococcaceae Phocea massiliensis 0.055
bin_Denmark.72 Bacillota_C Negativicutes Selenomonadales Selenomonadaceae Megamonas funiformis 0.082
bin_Denmark.85 Bacillota_C Negativicutes Selenomonadales Selenomonadaceae NA 0.077
bin_Brazil.5 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella sp900540415 0.080
bin_Malaysia.18 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotellamassilia sp000437675 0.051
bin_Brazil.48 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Bacteroides sp900766005 0.050
bin_Malaysia.117 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotellamassilia sp900541335 0.047
bin_CaboVerde.18 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella copri 0.038
bin_Aruba.10 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella sp900544825 0.037
bin_Malaysia.64 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Paraprevotella clara 0.033
bin_Brazil.163 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella lascolaii 0.029
bin_Spain.4 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola sp000436795 0.027
bin_Spain.48 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Bacteroides fragilis 0.019
bin_Brazil.96 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Bacteroides uniformis 0.014
bin_Malaysia.105 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola massiliensis 0.014
bin_Spain.21 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola coprophilus 0.012
bin_Malaysia.121 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Bacteroides stercoris 0.011
bin_Brazil.43 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola vulgatus 0.010
bin_Denmark.30 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola sp900546645 0.010
bin_Brazil.103 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola plebeius_A 0.009
bin_Malaysia.77 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola sp900542985 0.006
bin_Brazil.6 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola coprocola 0.005
bin_Brazil.110 Bacteroidota Bacteroidia Bacteroidales Barnesiellaceae Barnesiella intestinihominis 0.024
bin_Malaysia.13 Bacteroidota Bacteroidia Bacteroidales Muribaculaceae NA 0.075
bin_Spain.23 Bacteroidota Bacteroidia Bacteroidales Muribaculaceae CAG-279 sp900541935 0.051
bin_CaboVerde.37 Bacteroidota Bacteroidia Bacteroidales Porphyromonadaceae NA 0.058
bin_Brazil.11 Bacteroidota Bacteroidia Bacteroidales Rikenellaceae Alistipes putredinis 0.045
bin_Malaysia.131 Bacteroidota Bacteroidia Bacteroidales Rikenellaceae Alistipes putredinis 0.043
bin_Brazil.38 Bacteroidota Bacteroidia Bacteroidales Rikenellaceae Alistipes communis 0.039
bin_Brazil.86 Bacteroidota Bacteroidia Bacteroidales Rikenellaceae Alistipes dispar 0.030
bin_Brazil.138 Bacteroidota Bacteroidia Bacteroidales Tannerellaceae Parabacteroides sp000436495 0.023
bin_Brazil.124 Bacteroidota Bacteroidia Bacteroidales Tannerellaceae NA 0.019
bin_Brazil.160 Bacteroidota Bacteroidia Bacteroidales Tannerellaceae Parabacteroides distasonis 0.013
bin_CaboVerde.10 Campylobacterota Campylobacteria Campylobacterales Campylobacteraceae NA 0.070
bin_Brazil.68 Campylobacterota Campylobacteria Campylobacterales Campylobacteraceae Campylobacter_D helveticus 0.049
bin_Spain.26 Campylobacterota Campylobacteria Campylobacterales Campylobacteraceae Campylobacter_D upsaliensis 0.047
bin_Brazil.145 Cyanobacteria Vampirovibrionia Gastranaerophilales Gastranaerophilaceae NA 0.044
bin_Brazil.140 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_A sp900543175 0.029
bin_Malaysia.6 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_B sp900545035 0.015
bin_Brazil.159 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_B sp900541465 0.012
bin_Brazil.146 Pseudomonadota Alphaproteobacteria RF32 CAG-239 CAG-495 sp001917125 0.074
bin_Malaysia.50 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae NA 0.084
bin_Brazil.82 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum succiniciproducens 0.046
bin_CaboVerde.33 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum sp900543125 0.045
bin_CaboVerde.55 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum_A thomasii 0.036
bin_Brazil.111 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Succinivibrio sp000431835 0.024

8.3.2.2 Estimate - support plot

estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="sexMale") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
    dplyr::select(genome,mean)

support <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="sexMale") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
    dplyr::select(genome,support)

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
  filter(genome %in% estimate$genome) %>%
  arrange(match(genome, estimate$genome)) %>%
  dplyr::select(phylum, colors) %>%
  unique() %>%
  arrange(phylum) %>%
  dplyr::select(colors) %>%
  pull()

inner_join(estimate,support,by=join_by(genome==genome)) %>%
    mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
    mutate(support=ifelse(mean<0,1-support,support)) %>%
    left_join(genome_metadata, by = join_by(genome == genome)) %>%
    mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
    ggplot(aes(x=mean,y=support,color=phylum))+
      geom_point(alpha=0.7, shape=16, size=3)+
      scale_color_manual(values = phylum_colors) +
      geom_vline(xintercept = 0) +
      xlim(c(-0.4,0.4)) +
      labs(x = "Beta regression coefficient", y = "Posterior probability") +
      theme_minimal()+
       theme(legend.position = "none")

8.3.2.3 Predict responses

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = c("Male","Female","Unknown")
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "sex", 
                      non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
            predict(m, Gradient = ., expected = TRUE) %>%
            do.call(rbind,.) %>%
            as.data.frame() %>%
            mutate(sex=rep(gradient,1000)) %>%
            pivot_longer(!sex,names_to = "genome", values_to = "value")
# weights:  4 (3 variable)
initial  value 63.769541 
final  value 61.728471 
converged

8.3.2.4 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- pred %>%
  group_by(sex, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "sex") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="sex")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "sex") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)
8.3.2.4.1 Positive associated functions at element level
# Positively associated

unique_funct_db<- GIFT_db[c(2,4,5,6)] %>% 
  distinct(Code_element, .keep_all = TRUE)


element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  filter(positive_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)
# A tibble: 30 × 9
   trait    mean       p1      p9 positive_support negative_support Domain       Function                           Element       
   <chr>   <dbl>    <dbl>   <dbl>            <dbl>            <dbl> <chr>        <chr>                              <chr>         
 1 B0220 0.00893 0.00137  0.0165             0.924            0.076 Biosynthesis Amino acid biosynthesis            Beta-alanine  
 2 B0212 0.0293  0.00379  0.0522             0.922            0.078 Biosynthesis Amino acid biosynthesis            Arginine      
 3 B0218 0.0110  0.000319 0.0224             0.902            0.098 Biosynthesis Amino acid biosynthesis            Tyrosine      
 4 B0310 0.0183  0.00536  0.0312             0.964            0.036 Biosynthesis Amino acid derivative biosynthesis Tryptamine    
 5 B0303 0.0145  0.00204  0.0255             0.924            0.076 Biosynthesis Amino acid derivative biosynthesis Ectoine       
 6 B0307 0.0418  0.00475  0.0711             0.914            0.086 Biosynthesis Amino acid derivative biosynthesis Spermidine    
 7 B0901 0.00118 0.000115 0.00255            0.917            0.083 Biosynthesis Metallophore biosynthesis          Staphyloferrin
 8 B0105 0.0262  0.0106   0.0428             0.948            0.052 Biosynthesis Nucleic acid biosynthesis          ADP/ATP       
 9 B0104 0.0351  0.0102   0.0581             0.942            0.058 Biosynthesis Nucleic acid biosynthesis          CDP/CTP       
10 B0106 0.0150  0.00305  0.0267             0.921            0.079 Biosynthesis Nucleic acid biosynthesis          GDP/GTP       
# ℹ 20 more rows
8.3.2.4.2 Negatively associated functions at element level
element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)
# A tibble: 22 × 9
   trait     mean       p1        p9 positive_support negative_support Domain       Function                   Element           
   <chr>    <dbl>    <dbl>     <dbl>            <dbl>            <dbl> <chr>        <chr>                      <chr>             
 1 B0216 -0.0486  -0.0729  -0.0249              0.026            0.974 Biosynthesis Amino acid biosynthesis    Tryptophan        
 2 B0215 -0.0591  -0.0903  -0.0279              0.03             0.97  Biosynthesis Amino acid biosynthesis    Histidine         
 3 B0208 -0.0292  -0.0541  -0.00505             0.069            0.931 Biosynthesis Amino acid biosynthesis    Valine            
 4 B0209 -0.0283  -0.0555  -0.00322             0.073            0.927 Biosynthesis Amino acid biosynthesis    Isoleucine        
 5 B1012 -0.00833 -0.0130  -0.00368             0.005            0.995 Biosynthesis Antibiotic biosynthesis    Fosfomycin        
 6 B1004 -0.00460 -0.00875 -0.000640            0.081            0.919 Biosynthesis Antibiotic biosynthesis    Bacilysin         
 7 B0604 -0.0259  -0.0515  -0.000348            0.098            0.902 Biosynthesis Organic anion biosynthesis L-lactate         
 8 B0704 -0.0583  -0.0890  -0.0273              0.013            0.987 Biosynthesis Vitamin biosynthesis       Pantothenate (B5) 
 9 B0711 -0.0330  -0.0522  -0.0149              0.041            0.959 Biosynthesis Vitamin biosynthesis       Menaquinone (K2)  
10 B0710 -0.0179  -0.0297  -0.00596             0.057            0.943 Biosynthesis Vitamin biosynthesis       Phylloquinone (K1)
# ℹ 12 more rows
positive <- element_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- element_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
#      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
#      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.3.2.5 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- pred %>%
  group_by(sex, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "sex") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="sex")
   })
#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "sex") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
function_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()
tinytable_f111aufjkf4li5ao75t9
trait mean p1 p9 positive_support negative_support
D07 0.0375944479 1.972014e-02 0.056148600 0.992 0.008
B01 0.0120235113 2.924545e-03 0.021231683 0.939 0.061
B03 0.0134706279 3.129619e-03 0.021856061 0.932 0.068
D06 0.0046369186 1.904574e-04 0.009455340 0.909 0.091
B09 0.0008808116 2.773819e-05 0.001863172 0.905 0.095
D03 0.0086358818 -2.117773e-03 0.018844548 0.874 0.126
B04 0.0086819802 -1.265781e-03 0.019197647 0.866 0.134
D05 0.0037777625 -6.379054e-03 0.011495393 0.822 0.178
B06 0.0095293505 -4.760244e-03 0.023419874 0.812 0.188
D09 0.0044482653 -4.066690e-03 0.012324899 0.808 0.192
S03 0.0154206703 -1.668538e-02 0.041133011 0.800 0.200
D08 0.0003742537 -4.028928e-04 0.001116957 0.691 0.309
D01 0.0007972038 -4.891676e-03 0.006278674 0.668 0.332
# Negatively associated
function_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()
tinytable_5wmwplun1stiulfbo119
trait mean p1 p9 positive_support negative_support
B10 -0.0006304998 -0.0009540545 -0.0003081384 0.020 0.980
D02 -0.0151032927 -0.0273709242 -0.0026929023 0.061 0.939
S01 -0.0134498482 -0.0281931644 0.0021095325 0.132 0.868
B07 -0.0093150264 -0.0210303996 0.0025766178 0.142 0.858
S02 -0.0039566793 -0.0121727363 0.0061940272 0.170 0.830
B02 -0.0035889597 -0.0161455535 0.0073469073 0.362 0.638
B08 -0.0017867040 -0.0149536701 0.0101457248 0.533 0.467
positive <- function_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- function_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
#      xlim(c(-0.02,0.02)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")